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Abstract 


In this paper, a new orthogonal system of nonpolynomial basis functions is 
introduced and used to solve a class of time-fractional partial differential 
equations that have nonsmooth solutions. In fact, unlike polynomial bases, 
such basis functions have singularity and are constructed with a fractional 
variable change on Hahn polynomials. This feature leads to obtaining more 
accurate spectral approximations than polynomial bases. The introduced 
method is a spectral method that uses the operational matrix of fractional 
order integral of fractional-order shifted Hahn functions and finally converts 
the equation into a matrix equation system. In the introduced method, no 
collocation method has been used, and initial and boundary conditions 
are applied during the execution of the method. Error and convergence 
analysis of the numerical method has been investigated in a Sobolev space. 
Finally, some numerical experiments are considered in the form of tables 
and figures to demonstrate the accuracy and capability of the proposed 
method. 
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1 Introduction 


In recent works, science and engineering researchers found that the use of 
fractional calculus in modeling gives a more realistic description of various 
complex phenomena with long-range temporal cumulative memory. Frac- 
tional order operators have nonlocal and memory features. Therefore, these 
two important properties simulate and describe a variety of engineering and 
scientific problems with memory characteristics and inheritance more appro- 
priately than integer order differential equations, such as finance [30], physics 
[36], and hydrology [3], by using fractional differential equations. Many an- 
alytical methods have been used to solve fractional differential equations, 
such as the Green function method, Fourier, Laplace, and Mellin transform 
methods [24]. The complexity of integral and fractional differential operators 
and also the nonobservance of many properties expected in classic calculus 
encouraged researchers to study effective and reliable numerical methods for 
solving fractional differential equations. These numerical methods mainly 
include finite element and finite difference methods, spectral methods, and 
so on [31, 9, 8, 29, 34, 11, 1, 13, 22, 5, 23, 19]. In solving fractional or- 
der differential equations, two basic features that make classical methods 
not efficient and accurate are that fractional order operators have nonlocal 
properties and the other is the singularity of the solutions of fractional equa- 
tions. Therefore, spectral methods based on ordinary polynomials, which 
have high accuracy for solving problems with smooth solutions (see, for ex- 
ample, [33, 12]), are not suitable for solving fractional differential equations 
with nonsmooth solutions since they do not have high expected accuracy. 
Concerning the numerical solution of partial differential equations dependent 
on time, one of the most common approaches is to use the finite difference 
approximation together with the spectral approximation for time and spatial 
derivatives, respectively. One of the main drawbacks of this approach is that 
the temporal discretization error may overcome the spatial discretization er- 
ror, and the unknowns have to be solved simultaneously at all times [20]. As 
emphasized above, fractional differential equations mostly have nonsmooth 
solutions. It is also possible to encounter coefficients in terms of the given 
fractional equation in a nonsmooth case. On the other hand, in most of 
the spectral-introduced solving methods, in order to achieve high accuracy, 
they raise unrealistic assumptions. For instance, one of the assumptions in 
most of them is the smoothness of the unique regular solution of the frac- 
tional differential equation at the initial time t = 0 [32, 21, 35, 16]. So far, 
very few works have been done to solve fractional differential and integral 
equations with nonsmooth solutions, numerically, some of which can be seen 
in [15, 25, 26]. Due to their high accuracy, spectral methods have become 
one of the first choices researchers study to solve fractional differential equa- 
tions with nonsmooth solutions. Among these techniques, we can refer to 
the methods available in [37, 38, 7]. Analytical and numerical studies indi- 
cate the exponential convergence of these methods for nonsmooth solutions 
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in certain situations, and by using specific techniques, though, the exact so- 
lution of fractional time differential equations does not generally follow the 
mentioned form [27, 14]. 


This paper is organized as follows: fractional Hahn functions and their 
properties are defined in section 2, and also, function approximation and the 
operational matrix of fractional integration are introduced. In section 3, our 
proposed method is described. An error analysis is presented in section 4, 
and finally, some numerical examples are depicted in section 5. 


2 Fractional-order shifted Hahn functions approximation 


The main goal of this section is to introduce a new class of fractional basis 
functions, which are defined using shifted Hahn polynomials (SHPs) and 
applied to calculating their operational matrix of fractional integration. 


Definition 1. For given constants 01,02 > —1, and M € N, Hahn polyno- 
mials on [0, /] are defined as [17] 


~ ( k)i(k t O1 t 02 t 1)i( L)i 
hrp(x301,02,M) = ; , k=0,1,2,...,M, 
(x; 01,02, M) a (a, + 1);(—M),i! 
(1) 


where (-); is the Pochhammer notation, which is defined as 


(Co =1, (2) 
(Ci =C(C+1)---(¢+i-1), tEN, for¢eR*. 


Remark 1. The relationship between Stirling numbers and Pochhammer 
notation is as follows: 


(—h)s = (—1)h D7 527K, (3) 
1=0 
where 3) are Stirling numbers of the first kind defined as 


i-l 3 
(1) _4)" cL? 2r—1 (r) 
si d ) a i-l—r Fi—ltr) 


() 


in which s;’ are Stirling numbers of the second kind in the form 
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Now, by using (3) in (1) and the changing of variables as « = “4, we can 
achieve the following standard polynomial form of SHPs on (0, L] as 


_ Mt 
h(t; 01, 02, M, L) = he(—;01, 02, M) 


L 
ki 
1 M 
=> yi SS a x $0( yH4! 
é (o1 + 1);(—M)ji! £ 
i=0 1=0 
ki 
= Ai rat’, 
i=0 1=0 


for k = 0,1,2,...,M, where Aj. = (—1)'$ =e ol recess 


SHPs are orthogonal on [0, Z] via the inner product in the following form 
[28]: 


M 
(Fada = Fra Zr a(n), (4) 
r=0 


where w(r) is a real nonnegative weight function defined by 


ae — (ort 2\ fo. + M-« 
B(cior,oa,M) = ( . ( Vi ik (5) 
The orthogonal relationship of SHPs is as follows: 
M 7-2 ie : 
(hy hj) _ ye: hr, (Gr, 01,02, M, L)a(r), k =); (6) 
i: 0, k#& j. 


To define fractional-order shifted Hahn functions (FOSHFSs), t is substituted 
by t® in SHPs such that a is a positive real number. Therefore, FOSHFs can 
be defined in the following form: 


ass. )i(k +o, +02 +1); (1), M pat 
he (t; 01,02, M,L) = ou ee M);,i! | x 5; (Gt 
i=0 1=0 
ne k =0,1,2,...,M. 
1=0 1=0 


(7) 


Proposition 1. FOSHFs are orthogonal on [0, Z] via the inner product in 
the following form: 


=D (557) )al(5p7)* 20), () 
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where W(r) is defined in (5). 


Proof. Replacing f,g by he, nei in (8) and then using the orthogonality prop- 
erty (6), the assertion is available. 


Definition 2. Associated with the FOSHFs, the orthonormal FOSHFs can 
be defined as 


He (t; 01, 02,M, L) 
1 —_ 
= — h@(t;01,02,M,L). (9) 
/ PRE; 01,02, M, L), h&(t; 01,02, M, L))e 


2.1 Function approximation 


For an integer m > 0, the Sobolev space H” [a, b] is 
H®(a,b] = {u € L2 [a,b]: 0< 7 <m, u(x) € L2[a, b]}, 


where L?2 is the space of all square-integrable functions with respect to the 
weight function @. Indeed, H2"[a, b] is defined as the vector space of functions 
u € L2 [a,b] such that all derivatives of u of order up to m can be represented 
by functions in L2, [a,b]. 

Goertz and Offner described the expansion of a function by Hahn polynomi- 
als and concluded that the series expansion of a function by Hahn polynomi- 
als converges pointwise under some assumptions (for more details, see [10]). 
Therefore, any function u(t) € L2[0, L] can be expanded in terms of FOSHFs 
basis. In practice, only the first (IM +1) terms of FOSHFs are considered. 
Hence 


M 
t) = yy gH (t; 01,02, M, L) = um(t) = UTH™ (t; 01,092, M, D), (10) 
i=0 
where U" = [ug,uj,..., Uys] is the vector of FOSHFs coefficients, which can 


be derived as 


uw = me ), HE (t; 01,02, M, L))~ 


Be 


L D4 : 
= Soul ((a77) HE((sqr rye; 01,02,M, L) és (r ); i=0,1,...,44,(11) 


and H\%) (t; 01,02, M,L) is the vector of FOSHFs defined as follows: 


Ho) (ts 01,92, M, L) = [Hee 01,92, M, L), HE (t; 01,02, M, L), 


Iran. J. Numer. Anal. Optim., Vol. 13, No. 4, 2023,pp 672-694 


677 A shifted fractional-order Hahn functions Tau method ... 


., HS, (t; 01,02, M, L)]". (12) 


For simplicity, from now on, H(t; 01,02,.M, L) is presented by Ho (t). 
Similarly, any two variables function f(z,t) € L2([0,L] x [0,T7]) can be 
approximated by the FOSHEFs as follows: 


MN 
2) ae (2301, 02, N, L)H5 (t;01,02, N, L) 
7=0 7=0 
=: fu,w(2,t) = (HD (2) FHY (0), (13) 


where F' = [f;,;] is an (+1) x (N +1) matrix that its entries are 


T .as- Do va 
hs = (ri)*, (Gera) YAR (Sori) *s01,02,M, L) 
r1=0 r2=0 
#8 ((-rs)4 are 
Hy; ((sp72)? 301, 02, N, T)o(r1)o(r2), (14) 


for i=0,1,...,M and j =0,1,...,N. 


Theorem 1. Let M,N €N, A = (0, LZ] x [0,7] and let f € H2(A). Sup- 
pose that fi.n(z,t) = HO (0 ))? FU) (t) is the best approximation of 
finQ= span{HF (a; 01,02, M, LYH3(t;01,02, N,T)| 1=0,1,...,M, j = 
0,1,...,N}. We will have 


LMt+27M+2 


Fat) = Faw (oe, t)\lZ2 (a) S Q2M+N+1)(M +1)!(N +1)! B 


2s QMtN x,t) ae alt 
where F = max(z,t)ea |G | such that g(a, t) = f(@e,t@). 


Proof. Let dun (n,€) be the interpolation polynomial of g(n, €) = f(n=,€2) 
at (M + 1)(N +1) shifted Chebyshev points in A. Then 


Tywe max oer on §) 
2 (ngyea  OnMOEN © 


1 L 


Mang | 


la(n, €)—-oun(n, €)| < QM+N ( 


sb M+N 
If F = max(y¢yea |r| and n = «%, € = 1° are sets, then we get 


a a 1 LD T eS 
|g(x ie) — bu,n(x ,t?)| = Q2M+N(M a 1)! ( 2 yer 2 es (15) 


It is obvious that dn (x%,t®) € Q. So, since fiz. (z,t) is the best approx- 
imation of f concerning L? — norm, we have 


II f(a, t) — fat, (@, thlle < Ilf(@,t) — om,n (2, Alle 
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L eke 
= (f | (f(a, t) — ou.n (a, t)) dtdx)?. 


Thus, from (15) the assertion is derived. 


2.2 FOSHFs operational matrix of fractional integration 


Here, an operational matrix of fractional integration for FOSHFs is going to 
be obtained. Note that the Riemann—Liouville fractional integration of order 
6 for a function f is defined as 


1 x 
I’ f(x) = aa / (a — t)°-1 f(e)dt, z>a, v>0. (16) 
TO) Ja 
For this special type of fractional integration, there are some particular prop- 


erties. The most useful of which is 


y Toe} +y 
gy 0 
I" = Pept) r a ‘ (17) 


Using the above concepts, the following lemma states the FOSHFs oper- 
ational matrix of fractional integration. 


Lemma 1. The fractional integration of order 8 of the vector HH) (t) can 
be expanded by itself as follows: 


PH (t) ~ PoH (0), (18) 


where By = [pxj](ar+1)x (+1); Which is called the FOSHFs operational ma- 
trix of fractional integration with 


Mok GG i _ i 7 
Prj = »; ye > > Gr) Ave Aa aa hen (Gorse, 


r=0 i=0 [=0 113=01,;=0 


Proof. According to (12), we have 


I?H (t; 01, 02, M, L) 
IH? (t; 01,02, M,L) 


PHO (t) = (19) 


I°H?(t; 01, 02, M, L) 
PH (t;01, 02, M,L) 
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By using (7), Proposition 1, the linear property of operator I, and (17) for 
each entry in (19), we will have 


ki 
I°HR(t; 01,02, M, L) = SO Aare 


i=0 1=0 
T(al+1) os) 
=P YB pty pete 
i=0 I=0 (al +9 + 1) 


= Sa (t;01,02,M,L), k=0,1,...,M, 


where 
ki 
7a T'(al + 1) a Wa a 
Pj = oD Aik 'FaltetD Hs (t;01,02,M,L))3 (20) 
i=0 1=0 
M ko 4 
: or T(al + 1) L ee be 
= A; “((—r)e; M,L 
2H) DD, kl Mal+0+ mM" ry FG (aor r)* 301,02, ) 
Mik i 
= hn T(al +1) DL dla ——— DL at, 
— 2 O()ALnIT Ta 1+0+1) mM” H; (97 r) 301,02, M, L), 
r=0 i=0 1=0 
where Ajxy = Sie . substituting (7) in (20) 


VJ (he (t301,02,M,L),h2 (t:01,02,M,L))% 
instead of Ho ((4r)a ;01,02, M,L) finishes the proof. 


3 Description of method 


The main aim of this section is to approximate the solution of the following 
equation: 


D?u(a,t) = —au(a, t) +b an f(a, t), (21) 


subject to the initial and boundary conditions 


u(z,0)= g(x), u(0,t)=A(t), u(L,t)=n(t) forO<a¢<L,0<t<T. 


Let 
07 u(a, t) 
Ox? 
Applying the integration operator I” on both sides of (22) and using the 
operational matrix of integration (18), for 0) = 1 and J = 2, respectively, 
yield 


~ (H%(x)) UH (2). (22) 
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Ou) Hg, (0))? (BY UHR (t) + xh(0), (23) 


u(a,t) = (HGr(w))” B?) UH (t) + eh(t) + X(t), (24) 


in which the function h(t) is calculated by putting « = LD in (24) and then 
applying the final condition u(L,t) = n(t) as follows: 


h(t) = F(nlt) — A) — S4(L))" (BYP UHR (1). 


Therefore (25) and (24) can be rewritten as follows: 


Pu) Cage a) OB UHR() 
+ + (nlt) — A(t) — (HE_(L))? (PB?) UKE), (25) 

u(x,t) ~ (H9,(x))? (PB?)T UMA (t) 
7 (nt) — A(t) — (Hi (L))" OR?)TUHK(O) +A). (26) 


By substituting (22), (25), and (26) in (21), we get 
Deu(a, t) ~ —a[(HR;(a))? OB?) UH (t) 
+ F(n(t) — A(t) — (HRr(L))" 9B?) UA () + 0) 
+ b[(HGr(a))? (PB?) UNA (t) 
~ = (nlt) — A(t) = (HR,(L))" (B?)" UH (t))] 
+ c(HGp(a)) UH (t) + f(@, t) 
= (H§,(x)) AHA (t), (27) 


where 
= —a(P?)TU + 0% (HG (LE)? R)TU 
+ DOB)PU — b= (HGi(Z))" B)U + U + Ka, 


and K,, X, and 1 are the matrix and vector coefficient of FOSHF-approximation 
related to the following relations: 


ky (x,t) = f(a, t) a(F(n(t) d(t)) + A(t)) + = (n(t) — A@)) 
~ (HGr(x))"KiHy (t), 


w ~ (Hf;(x))"X, 
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T 
Le (Hp(a))"1. 


Again, by applying the integration operator J? on both sides of (27) and 
using the operational matrix of integration 8”, we will have 


u(a,t) ~ (Hfp(a)) WP" Hy (t) + g(2). (28) 


Equalizing the right sides of (26) and (28), we get 
x 
(Hiir(2))" (0B?) UF (Hiir(L))" (B?)* U-258"}Hy (4) = (Hf; (x))" Koy (t), 


where Kg is the matrix coefficient of FOSHF-approximation related to the 
following relation: 


ka(2,t) = g(x) — F(n(t) — A(t) — AW) = (HG (@)" KoA (0. 


Thus x 
?)7U— 5 (Hiie(L)) (BU — AP" (t) = Ko, 


which can be rewritten as 
BU + CUD = E, (29) 


where B = (I — F(H§,(L))")(P?)*, 

€ = [al — 2X(HG,(L))" + 21(H9,(L))7 12)" — WPT — cl, D = P", and 
€ = K;$" + Ky. Equation (29) is a matrix equation with the unknown 
matrix U. It can be solved by the global GMRES method. After solving the 
equation, by placing the obtained matrix U in (24), the approximate solution 
of the problem is obtained. 


4 Error analysis 


In this section, the convergence of the introduced method in a Sobolev space is 
considered. An upper bound is derived for the absolute error of the proposed 
method. To this end, some bounds are obtained for the approximations of 
different parts of the mentioned equation. First, the basic definitions and 
concepts related to Sobolev spaces are from the books [4, 18], with a slight 
change in symbols. 


Let A be an open subset of R” and let L2(A) be the space of all square- 
integrable functions concerning the weight function w. For an integer m > 0, 
the Sobolev space H2’(A) is 


Ae (A) = {u] we EA (A), OvYUE EA) for all |v| < m}, 
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where 0” is called the distributional derivatives and defined as the following 
form: 


avy 
= ? 
Oxy! 0x5? ++ - Oxn” 


For all u, v € H2°(A), the inner product is given as 


O”’u 


lary Hig oh ity 


(U,V) am (ay = (U,V) £2.(A) + S> ( (O"u, O"V) 12 (A): 


1<|v|<m 


The corresponding norm and seminorm are defined as 


lela cay = (ellZzcay + SO (O’ullzecy)?: 


ie 


lulae (ay = (20 \|o” ullz2 (i)? 


It is obvious to see that if m > 0, then |[ul]z2 (a) < ||ullaz(ay. In a special 
case, for m = 0, it yields ||ull~(a) = |lullz2 (a). Also, for m = 0, we have 
[ule (ay = |lullz2 cay- 

Suppose that u € H2"(A) and Poy are the orthogonal projection opera- 
tor, where A = [0, L] x [0,7] and 


M N 


Pate = yy ane (x; 301,02,N, LYH§ (t; 01, 02, N, L). 
1=0 j7=0 


In other words, Pri =um,n(2,t) = (Ho (x) TUHO (t). 
In the following, for simplicity and brevity, M = N, a = 6, and Py := 
Pow are stated. According to [6], for all u € H2(A), we have 


Ju — Paral za cay S CMP ™ ul rmiateny, OS FSM, 0) 


y? 
where C is a constant independent of M and only depends on m, 
’ 0, j =9, 
i= {oa j>0, 


and 


m 2 
7 i 
lly =( D> DDS) 2. 


k=min(m,M+1) i=1 


Theorem 2. Suppose that u(x,t) © H3’(A), m > 0 and that uy, ac(z, t) is 
the best approximation of u. Then 


Ilu(a, t) — una (2, 4)Il22 (ay S lu(@, 2) — ww ae (@ Mls ay (AD) 
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<COMPO—™ul pmiveay, OSG <m. (31) 


Proof. Considering this fact that || - |[72(a) < || - llama), imequality (30), and 
the uniqueness of the best approximation, the proof of the theorem is easily 
done. 


Lemma 2. Suppose that the assumptions of Theorem 2 are true, that 
u(x,t) ~ uynv(2,t) = (HO (x))TUHO(t), and that By is the FOSHF- 
operational matrix of fractional integration. Then 
| ou(a,t) — (Hyp (@)) BSUHY Olay 
L? 


es p(j)—m : 
<To+pe [el mt (ays 0<j<m. 


Proof. According to (16), we have 


lZ2u(w, t) — (HY? (w) )PRSUHY? (#)lhzz ca) 
= |[IPu(z, t) — umn lacy 


= ||Ip (u(x,t) — umn (2, t))lIr2 (a) 


='|| rw) [e — €)°"*(ul€,t) — umn (E,t))d8 || 2 (a) 
= ral * (u(a, €) _ um,n(2,t))||12(a)- (32) 


Now, by using this fact that ||f * g\|p < ||fIl1-llg||p, and Theorem 2, respec- 
tively, we get 

[® 462) Na T yyy) (t)|| 

|Z u( ( M (x))" Py N L2 (A) 


x,t) — 
L? 
< orca (4) — um (2, t)|| 52 (A) 
L® 


IA 


(j)—m ; 
ro - Fl J [ul mime (ays 0 < J <m. (33) 


To get an error bound for derived approximation in the proposed method, 
which has been introduced in section 3, without losing the generality, we 
suppose that 


a > (HP (x) TUH® (t) = oun (2, ¢), (34) 
oulest) ~ (HY (x) WHY (t) = oun (2,t), (35) 
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u(a,t) ~ (HO (x) VHY (t) =: vw (2, t)- (36) 


As it can be seen from the process of the presented method in section 3 that 
relation (26) has appeared in applying the operator I? on the sides of (34) 
and (36) can be derived by expanding (26) in terms of FOSHFs basis. It is 
easy to see that 


0 u(z, t) 


—_ 2 
I|u(x, t) — bun (2, t)[l22 (ay = 2 —a a 


— Foun (a, #)|l 22 (ay: (37) 
So, considering (37) and applying Lemma 2, the following corollary is ob- 
tained. 


Corollary 1. If relation (34) is true, then 


07 u(z, t) 


Aut le (ays 0 < j <m. 


i? 4 aps 
I|u(z, t)—bu,n (2, tI] 22 (ay < maou | 


Consider the main equation (21) and the presented method in section 3, 
by substituting (34)—(36) on the right side of (21) and applying the operator 
I? on it, we get 


u(x,t) ~ —al/’h(a,t) + bly pa, t) + ely O(a, t) + I? f(x,t) + 9(z). (39) 
On the other hand, we have 


2 
1) =r ea a lr t) 


on tT? f(x,t) + g(a). (40) 


Putting the right side of (39) and (40) as equivalent, we define perturbation 
term as follows: 


Ou(x,t) 
Ox 


— oun (a,t)). (41) 


Run (a, t) = —al? (u(a,t) — bun (a, t)) + v1? ( 


07u(z, t) 


— yu,n(a,t)) + cl} ( a2 


Theorem 3. Suppose that, u(x,t) € H2’(A) for m > 0 is the exact solution 
of (21). If dyr,n(a,t) is the approximate solution, obtained by applying the 
presented method, then Ryn (x,t) —+ 0 as M, N — oo. 


Proof. According to (41), we have 


Razn (a, t)Ilzacay S lallly (u(@,t) — daw (2, t)) Ila cay 


Ou(x,t 
ee ee — ong wes t)) lagen 


OP u(a, t 
+ lene EE — baaw(ePllugey. — 2) 
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Now, by applying Lemma 2 in approximations (34)—(36), respectively, we get 


|Z? (u(a, i — om,n(2,t))||22(A) 


ats Ta OMe lulamatcys OS FSM, (43) 
ex one 2 bya (ast)loaca 
= TF tray wy Age; (44) 
ee eerteiilaw 
< saenomn ie may OSJSm. (45) 


So, by using (43)—(45) in (42), it yields 


|v (2, t)|] 22 (a) 
re oe du Ou 
< Tornoe (la|fetl prmsae cay + [Pll a re cay ae lell alm (ay): 


Hence, it is concluded that Ryz,n (x,t) —+ 0as M, N —> co. Oo 


5 Numerical results 


In this section, the introduced method in section 3 is utilized to approximate 
the solutions to problems. It should be mentioned that the maximum of 
absolute error is the infinity norm of the error function and 


Lio = ;, L’)|. 
so = max, le( 2,7) 


All numerical experiments have been performed using MATLAB R2017a on 
a Core(TM)2 laptop with 4GB RAM and a speed of 2.00 GHz. 


Example 1. Consider the following time-fractional diffusion differential 
equation: 


u(x, t) 


Doula, t) = —u(a,t) + F2 


+f(z,t), O0O<0<1, (2,t) € [0,1] x [0,1], 

(47) 

where f(x,t) = sin(rax)(1 4 wi wee. subject to the initial and 
boundary conditions: 


u(z,0)=0, u(0,t)=0, u(1,t) =0. 
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The exact solution is u(x,t) = roan sin(r2). Table 1 shows the L,,-norm 
of absolute error for fixed N = 1 and some M and £ in comparison to [2]. In 
Figure 1, the L..-norm of absolute error for fixed N = 1, J) = 0.9, and some 
M =4,5,...,10is shown, which demonstrates that the approximate solution 
converges to the exact solution as M increases. Finally, Figure 2 shows the 
absolute error functions for fixed N = 1, J) = 0.9, and some M = 6, 8, 10. 


107 


105 


L,,-Norm of Absolute Error 


10% 


10° 


4070 


Figure 1: Loo-norm of the absolute error function for fixed N = 1, J = 0.9, and some 
M =4,5,...,10 (Example 1) 


M=6 M=8 M=10 


Absolute Error 
Absolute Error 
Absolute Error 


Figure 2: Absolute error functions for fixed N = 1, 8 = 0.9, and M = 6, 8, 10 (Example 


1) 


Example 2. Consider the following inhomogeneous fractional-order Burger’s 
equation: 
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Table 1: Loo-norm of absolute error for fixed N = 1 and some M and # in comparison 
to [2] (Example 1) 


M 9 =0.25 [2] 9 =0.5 [2] 9 = 0.75 [2] 
4 2.301e-3 1.690e-3 2.389e-4 4.979e-3 2.287e-4 2.918e-3 
6 3.638e-5 5.764e-4 3.721le-6 3.33le-5  3.589e-6 2.75265 
8 4.517e-7 1.761e-6 4.621e-8 1.754e-7 4.455e-8 1798e-7 
10 7.101e-10 3.127e-9 7.263e-10 8.428e-10 7.003e-10 8.116e-10 
2 
0 O*u(az,t) Ou(za,t) 
Hieg - f(a,t), 0<9<1, (a,t) € (0,1) x [0,1], 
Ox Ox 
(48) 
2-0 « Bic fag tes “ype 
where f(a, t) = Facey + 22-2, subject to the initial and boundary conditions: 


u(x,0)=a27, u(0,t)=t, u(1,t)=14+#. 


The exact solution is u(x,t) = x? + t?. Figures 3 and 4 show the absolute 
error functions after solving the problem by using the presented method with 
M=2, N=4, a=1, @ = 0.5 for the fractional-order derivative 0 = 0.5 
and M=2, N=2, a=1, 8 =1, and J = 1, respectively. 


Absolute Error 


Figure 3: Absolute error function for M = 2, N = 4, a 1, 8 = 0.5, and 0 = 0.5 
(Example 2) 
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Absolute Error 


Figure 4: Absolute error function for M = 2, N=2, a=1, 8=1, and 0) = 1 (Example 
2) 


Example 3. Consider the following transformed time-fractional Black— 
Scholes model with homogeneous boundary conditions: 


2 92 2 
4 a* Ou(a,t) a Ou(x, t) _ 
D; u(x, t) ) Ox2 (r 2 ) Ox =r ru(x, t) ~~ fiz, t), 
0<0<\1, (2,t) € (0,1) x (0,1), (49) 
where 


(50 — 4e*) — r(a® — 2*)), 


subject to the initial and boundary conditions: 
u(x,0)=a2°—2*t, u(0,t)=0, u(l,t)=0. 


The exact solution is u(a,t) = (t? +1)(x° — 2+). Let r = 0.02 and let o = 0.8. 
Figure 5 shows the absolute error function obtained by applying the presented 
method for 0 =0.5a=1, 6 =0.5, M=5, and N =6. Also, Figure 6 shows 
the absolute error after solving the problem by using the presented method 
with M=5, N=83, a=1,and 6 =1 for v=1. 
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Absolute Error 


Figure 5: Absolute error function for ) = 0.5 a= 1, 6 = 0.5, M = 5, and N = 6 
(Example 3) 


Absolute Error 


Figure 6: Absolute error function for ? =0.5a=1, 8=1, M=5, and N =3 (Example 
3); 


Example 4. Consider the following time-fractional equation: 


2 
Dpu(x,t) = —au(e,t) + 22D 4 PMD _ F¢~0, 
0<¥<1, (a,t) €(0,L) x (0,T7], (50) 


subject to the initial and boundary conditions: 


u(a#,0) = x, u(0,t)=0, u(L,t) = Lee 
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where in the case of 0 = 1 and the function f is chosen as f(z,t) = 
—e'3 (x2 + 3\/x), the exact solution is u(z,t) = e~'a?. It is notable that 
in other cases of 0 < J < 1, the exact solution is unknown. Figure 7 shows 
the absolute error functions obtained by applying the presented method for 
0=0.5, a=0.5, B=1, M=5, and N = 4, 6, 8. Also, Figure 8 shows 
the Lo.-norm of the absolute error function for fixed M = 5, J = 0.5, and 
some N = 2,3,...,8, which demonstrates that the D..-norm of the absolute 
error function converges to zero as N increases. Finally, Figure 9 depicts 
approximate solutions for different 0< 0 <1, M=5, N =7, which shows 
that as 0 + 1, the approximate solution converges to the exact solution when 
v=1. 


Noa N=6 N=8 


10% x10® xo 


Absolute Error 
Absolute Error 
Absolute Error 


Figure 7: Absolute error functions for 0 = 0.5, a=0.5, 8=1, M=5, and N =4, 6, 8. 
(Example 4) 


L_, Norm of Absolute Error 


rs 


Figure 8: L.o-norm of absolute error function for fixed M = 5, 0 = 0.5, and N = 
2,3,...,8. (Example 4) 
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0.018 


0.016 


0.012 


0.008 


0.006 L L L L L L L L L J 
0 


Figure 9: Approximate solutions for different 0<0<1, M=5, N=7. (Example 4) 


6 Conclusion 


In this paper, a new orthogonal system of nonpolynomial basis functions, 
named FOSHFs, has been introduced and used to solve a class of fractional- 
time partial differential equations with nonsmooth solutions. For introducing 
the method, an operational matrix of fractional order integral of Hahn func- 
tions has been used for the first time as basis functions here. Furthermore, 
the convergence of FOSHFs approximation has been proved. In numerical 
examples, the obtained results have demonstrated the efficiency and conver- 
gence of the proposed method for the cases of nonsmooth solutions. 
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